Shearing of loose granular materials: A statistical mesoscopic model 
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A two-dimensional lattice model for the formation and evolution of shear bands in granular 
media is proposed. Each lattice site is assigned a random variable which reflects the local density. 
At every time step, the strain is localized along a single shear-band which is a spanning path on 
the lattice chosen through an extremum condition. The dynamics consists of randomly changing 
the 'density' of the sites only along the shear band, and then repeating the procedure of locating 
the extremal path and changing it. Starting from an initially uncorrelated density field, it is found 
that this dynamics leads to a slow compaction along with a non-trivial patterning of the system, 
with high density regions forming which shelter long-lived low-density valleys. Further, as a result 
of these large density fluctuations, the shear band which was initially equally likely to be found 
anywhere on the lattice, gets progressively trapped for longer and longer periods of time. This 
state is however meta-stable, and the system continues to evolve slowly in a manner reminiscent of 
glassy dynamics. Several quantities have been studied numerically which support this picture and 
elucidate the unusual system-size effects at play. 



I. INTRODUCTION 

Modeling the rhcology of granular media using contin- 
uum solid mechanics, has reached a high degree of sophis- 
tication in terms of constitutive equations . Whatever 
the complexity of load paths being studied, an accurate 
account of the experimental stress strain relationship can 
now be achieved provided enough parameters or internal 
variables are included in the constitutive laws. However, 
such approaches are descriptive and leave unanswered 
questions pertaining to the scale of grain sizes. 

In parallel to such phenomenological descriptive theo- 
ries, a lot of effort has been spent in recent years in devel- 
oping powerful computer models able to simulate granu- 
lar systems at the individual grain level. Molecular dy- 
namics approaches B , or other techniques such as "con- 
tact dynamics" || ||5 now offer the possibility of dealing 
with several thousand particles, and provide extremely 
realistic pictures of the detailed micro-mechanics. 

Such numerical techniques can be used to accurately 
investigate displacement fields, resolved both spatially 
and temporally 0]. The latter reveal an intriguing fea- 
ture: namely that even in the most simple tests, such as 
a simple steady shear imposed over large strains, the local 
displacement appears very unsteady, with short quiescent 
periods where the displacement field is spatially smooth, 
separated by sudden changes where the configuration of 
grains reaches a local instability and undergoes a rapid 
reorganization through significant displacements at the 
grain level. This temporal variability manifests itself in 
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giant stress fluctuations observed experimentally when 
particles and walls are stiff, and when the high frequency 
part of the stress signal is not filtered out||. Numerical 
simulations indicate |?J that instantaneous strain fields 
consist essentially in localized strains occurring along one 
or a few shear bands. However, as the strain increases 
(over the moderate range accessible in the simulation), 
there seems to be little or no correlation between suc- 
cessive shear bands, so that the time average of the dis- 
placement field erases these discontinuities and produces 
smooth strain fields. 

Such fluctuations are obviously ignored in continuum 
modeling. And indeed it may appear that the identifica- 
tion of these instabilities is relevant only for discussing 
fine details of microscopic and transient features. How- 
ever, their relevance can be judged only at a mesoscopic 
level of modeling, since the microscopic numerical tech- 
niques are far from being able to reach the relevant time 
scales. In fact, in the following we will argue that these 
instabilities may have a significant impact, both on large 
scale heterogeneities of the medium itself, and on a sys- 
tematic slow time evolution of the macroscopic friction 
angle. A short account of some of our results has ap- 
peared in a previous publication ||. 

The paper is organized as follows: in Section II, we 
will recall some features observed experimentally or nu- 
merically that we consider essential, and, in Section III, 
we progressively introduce the rules of a model whose 
aim is to describe some statistical aspects of shearing of 
loose granular media over large strains. In Section IV, 
we present in details the different quantities studied nu- 
merically for this model. We conclude in section V with 
a summary of our results and a discussion on possible 
experimental checks. 
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FIG. 1: Schematic picture of the shear process. The shear 
band is parallel to the shear direction z due to the periodic 
boundary conditions in this direction. We sum up along this 
direction to get a two dimensional sample in the xy plane. 



II. THE SHEAR PROCESS IN LOOSE 
GRANULAR MATERIAL 

We will address here the question of the behavior of 
granular media subjected to a simple shear for large 
strains. We restrict ourself to the simplest granular 
medium one may consider, namely rigid (undeformable) 
grains with Coulomb friction. This refers experimentally 
to dry sand subjected to a low confining pressure. We 
are concerned here with large strains, and thus in or- 
der to avoid the problem of boundary conditions which 
would limit the maximum strain, we consider an annu- 
lar shear cell. To simplify the problem further, we con- 
sider only the case where the problem is invariant along 
the shear direction. As shown in Figure |], the displace- 
ment is a single function of the coordinate of a radial 
cross-section (x, y), and constant along the orthoradial 
direction z (traditionally this situation is termed "anti- 
plane"). Moreover, we are interested only in the qua- 
sistatic regime, i.e., time as such is irrelevant, and only 
the total strain matters. Thus, in what follows, what is 
referred to as "time" is to be seen here as a practical 
means of parameterizing the total strain being imposed 
on the medium. 

One of the important observations of soil mechanics 
concerning such media is the concept of a critical state 
H 0. Depending on the preparation of the sample, 
the behavior under shear may differ considerably. For 
loose sand, (low density), the deviatoric stress to be ap- 
plied increases with the total shear strain and simulta- 
neously, a densification is observed jll], Q. However, as 
the shear strain increases, the density and shear stress 
seem to reach a plateau independent of the initial den- 
sity. This state is called the "critical state" . On the con- 
trary, if the initial density is large, a single shear band 
forms, while the rest of the medium remains frozen [fl3[ . 
The formation of the shear band is preceded by a vol- 
ume expansion of the medium |l4| , but after the band is 
formed, all further properties remain quasi constant. A 
detailed experimental investigation has revealed 15 that 
inside the shear band, the density tends to approach the 



critical state density. This concept of the critical state 
has received considerable experimental evidence over the 
years, and is implemented in a number of continuum con- 
stitutive laws. Experiments however mostly deal with a 
rather moderate total strain well below unity. 

A simple picture which is consistent with the critical 
state concept is that both the friction and the dilation an- 
gle increase with the density, and that the critical state is 
the density for which the dilation angle is zero (no change 
in volume under shear). Retaining the density as the only 
internal variable is an approximation. Other character- 
istics of the texture of the medium such as the fabric 
tensor (which has information about the orientation of 
contact normals), certainly play a significant role. For 
the purpose of simplicity, we will in the following only 
retain one single scalar internal variable governing the 
friction angle. It could be either simply the density or a 
combination of density and texture. Nevertheless, in all 
cases we will refer to this internal variable as "density" , 
irrespective of its precise meaning. 

As mentioned above, numerical simulations seem to re- 
veal evidence for the existence of instantaneous shear 
bands even in loose granular media. On the other hand, 
in experiments the strain appears to be homogeneous and 
not localized. The resolution of this apparent paradox is 
that the shear bands change rapidly, and may visit the 
entire medium in the process. Thus during an increment 
of shear, which can be observed experimentally, only a 
time average over many such shear bands is seen. In 
our modeling, we introduce a basic time scale for each 
elementary procedure. This time step is then clearly 
much shorter than most experimentally accessible time 
scales. However, our model attempts to achieve a qual- 
itative rather than a precise quantitative mapping. One 
of the main features of the model is to show that these 
two apparently unrelated facts: the existence of instan- 
taneous shear bands at early times and its localization 
at late times, are actually related, with a slow transition 
between these two limiting cases. This slow dynamics 
is reminiscent of slow ageing properties encountered in 
glassy systems, and indeed, we will see that a breakdown 
of ergodicity does appear in this model. 



III. THE MODEL 

A. Motivation and definition 

At every instant the two-dimensional medium is char- 
acterized by a single, scalar internal variable, the density 
g(x,y). This represents an average of the density along 
the orthoradial direction z. From this density, we deduce 
a corresponding local friction coefficient /i(i£, y). The lat- 
ter is assumed to be a single monotonically increasing 
function of the density |l(| . For simplicity, we may as- 
sume a linear relationship in the following although this 
is inessential. 

The strain is imposed on the shear cell through pre- 
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scribed displacements of the bottom and top planes. As 
particles are considered rigid, (no elastic deformation), 
the shear cell can only move if the shear force exceeds 
a threshold value proportional to the normal pressure. 
This limit stress is given by the "weakest internal sur- 
face" . Indeed, in our anti-plane geometry, the shear 
strain will localize on the surface (i.e. path in the (x, y) 
plane) which will fail first. The latter is assumed to be 
given by the following algorithm. For each directed path 
V spanning the entire cross-section along the a;-axis, we 
compute the maximum shear force it can support accord- 
ing to the local density. Assuming that the local slope 
of the path is always small, this maximum force F(V) 
is simply proportional to the sum of local friction coeffi- 
cients, and thus making use of the assumed linear vari- 
ation of the friction coefficient with the density, F(V) is 
proportional to the sum of local densities, 

S{V)= e(x,y), (1) 

(x,y)eV 

where the sum runs over the sites along the path. Among 
all the possible paths, the weakest V* (for which S(V*) 
min ; will fail first, and this fixes the value of the shear 
force F = F(V*). In agreement with the previously men- 
tioned observation, at every basic time step, the shear 
strain is realized along a single shear band. Away from 
this shear band, the strain rate is zero, and thus the den- 
sity is kept constant in time. However, inside the shear 
band, there is a relative motion of grains, and thus the 
density is susceptible to evolve. 

The next step is now to determine how the density in- 
side the shear band evolves with time. Even though it 
is observed that large strains are necessary in order for 
a system to reach its critical state, we argue that at the 
microscopic level, the evolution of the medium cannot 
depend on the total imposed strain. Thus the evolu- 
tion rules for the density within the shear band should 
be designed in such a way that they do not depend on 
the past history, but only on the present state (density 
field). As the density ought to contain the basic in- 
formation of the local characteristics, we propose that 
within the shear band, the local density g(x, y) is ran- 
domly modified. More precisely, in one elementary time 
step, corresponding to the "life-time" of the shear band 
in a very loose granular sample, we assume that the den- 
sity along the shear band acquires random uncorrelated 
values picked from a statistical distribution p(g). The 
uncorrelated character of the distribution is however jus- 
tified only on a mesoscopic scale. 

After the elementary strain event, we have a new den- 
sity map g(x, y). We now simply reiterate this procedure 
as long as desired: Namely we identify the new path 
which minimizes S(V), and update the value of the den- 
sities along this path randomly. As the purpose of the 
present article is to illustrate some statistical aspects of 
this dynamics, we do not try to mimic any specific gran- 
ular system by imposing a realistic density distribution 
or initial correlations in g. We will choose here a simple 




FIG. 2: a) Visualization of the tilted lattice. The shear band 
is marked with a thick line, b) A sample configuration of 
the minimal path on the normal square lattice with densities 
assigned to sites. 



uniform distribution between and 1 for p{g). The mean 
value and variance of the distribution p can be chosen ar- 
bitrarily, since a translation and rescaling of g does not 
affect the result. 

A key assumption of our model which may appear as 
precluding the occurrence of a slow evolution toward a 
critical state is the selection of the density values within 
a shear band from uncorrelated, smooth distributions. In 
fact, we will show below that, on the contrary, a collective 
and purely statistical effect produces a slow increase of 
the mean density over large strains. 

Our model is furthermore discretized on a regular 
square lattice. We have looked at two different kinds 
of square lattices to check the robustness of our results. 
In the first one the value of the density g is carried by 
the bonds. The orientation of the lattice is chosen so 
that the principal directions lies at ir/4 with respect to 
the (x,y) axis as shown in Fig. |] a). In the other ver- 
sion, density values are assigned to the sites of a square 
lattice. In this case the minimal path can be connected 
through the next nearest neighbours too as shown on Fig. 
U b) . Both square lattices give exactly the same results 
so in the following we just refer to them as square lattice 
realizations. 

We mention here that a third type of lattice, the hi- 
erarchical diamond lattice |l7j] was also studied. The 
numerical results are surprisingly, essentially unchanged 
by the unusual topology of this recursively constructed 
lattice. The easier construction of this lattice however 
allows us to solve the model analytically thus giving us 
a quantitative picture of the behaviour of the system. 
These results are presented in [fl8f where the intimate 
relationship of our model to other models of statistical 
physics is also discussed. 

The rules of our model, finding the extremal directed 
spanning path at every instant, is similar to finding the 
ground state of a directed polymer in a random potential 
|19| . However, in our case this potential is uncorrelated 
only at the beginning; it changes in time through the 
process described above, of ascribing new densities to all 
sites along the minimal path. It is clear from this re- 
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configurations and achieve a significant local densifica- 
tion. We also see within these wide and dense chan- 
nels, a single path with a lower density. This has been 
the last active minimal path in the channel. As time 
passes, the number of channels increases, and so does 
their width. They get partly interconnected, leaving al- 
ways the same scars of low density paths. Finally at the 
latest time shown on the figure, the average density is 
quite high, and traces of ancient minimal paths are still 
visible. Nevertheless what is striking is the occurrence of 
islands entrapped by these high density channels, where 
the density map looks like at the very early stage of the 
time evolution. This signals that these regions have ba- 
sically not been visited by the minimal path during the 
entire history of the system. These features are quite 
generic, and they reveal that the spatial (and temporal) 
organization of the activity is rather complex. The rest 
of the study is devoted to a more quantitative account 
of this activity, of the resulting kinetics of compaction, 
and the unexpected finite size effects which appear in this 
problem. 

In the following subsection, we will introduce the main 
measurements performed numerically on the model. 



FIG. 3: Snapshots of densities of a square lattice of size 256 by 
256 at time a) 1000 = AL, b) 10 4 = 40L, c) 10 5 = 400L and 
d) 10 6 = 4000L. The gray scale is presented at the bottom. 
The actual shear band is drawn in black. 



lationship between the models, that the shear band is 
sclf-afhnc with a Hurst exponent C, = 2/3 at the begin- 
ning, i.e. the transversal fluctuations of the band grow 
with the size (L x ) of the sample width as L|. We will 
discuss the time evolution of the roughness later in this 
paper. 



IV. NUMERICAL RESULTS 

We first show the density map, Fig. |], of the system 
at different times, t/L ranging from 4 to 4000. The grey 
scale chosen focuses on the vicinity of 1 so as to highlight 
the progressive densification. It may appear counter- 
intuitive at first that the rest of the medium shows a 
densification at all, when the only dynamics consists of 
finding a minimal path and updating sites along it ran- 
domly. However the reason is simply that this update 
systematically hunts out the sites with the lowest den- 
sity values and replaces them. 

At early times, we observe an apparently uncorrelated 
field. However, as time proceeds, it is possible to distin- 
guish preferential channels of high density aligned along 
the direction (x-axis) of the minimal paths. These chan- 
nels however have a significant width which shows that 
though the minimal path has been confined to this zone, 
it has enough freedom to explore different neighboring 



A. Definitions of numerically measured quantities 

1. Average density 

The most important quantity is the average density of 
the sample that we define as the mean of the density of 
the inactive sites, i.e. the sites not belonging to the shear 
band. We denote this by (g). This definition is conve- 
nient because (g) increases monotonically by the rules 
of our model. In experiments one of the most frequently 
measured quantities is the volumetric strain which is just 
the change in the inverse of the average density. In our 
model, when p(g) is chosen to be the uniform distribu- 
tion in the interval [0, 1], the average density is bounded 
by ((g) < 1). It is easy to see that in finite systems the 
steady state (the asymptotics) cannot be other than a 
system with maximal densities everywhere except for a 
path which will be always chosen as the minimal path. 
This state is equivalent to (g) = 1. 

We are interested in the approach to this asymptotics 
so we plot the quantity (1 — (g)) as a function of time. 
Within the granular medium context, this means that we 
mainly study a loose initial state and its convergence to 
the critical state |)| However, we will also present re- 
sults obtained when one starts from a high initial density 
later (Section |VG|). 

Figure |4| shows that the difference of the average den- 
sity from its asymptotic value first remains almost con- 
stant during a first stage t/L <C 1, and then it decreases 
steadily. This first increase of (g) is well captured by 
a reduced time equal to t/L. However, as time pro- 
gresses, the average density increases more and more 
slowly. Quite strikingly the larger the system size, the 



5 




t/L 

FIG. 4: The difference of the average density from its asymp- 
totic value 1 is plotted as the function of time. The system 
sizes are L = 32, 64, 128, 256, 512. The average was done 
over the inactive sites on the lattice and for an ensemble of 
20 to 1000 samples. 

slower the increase in density. Further down this will be 
interpreted as a breakdown of ergodicity. 



2. Shear band density 

It is also natural to define the density of the shear band 
that we denote by qsb- This is just the average den- 
sity of the sites along the minimal path (before updating 
them). As already mentioned in the introduction, we 
assume that the maximal static shear force is a single 
function of the density. Thus the density of the shear 
band can be related to the shear stress in experiments. 

Figure || shows the evolution in time of the difference 
(0.5 — qsb)- As expected, as time proceeds, the density 
along the shear band will tend toward the average of the 
random densities which are used to refresh the sites or 
bonds along the shear band. Using a uniform distribution 
of densities between and 1 implies that this average 
is 0.5. Note that in contrast to the previous case, the 
reduced time t/L accounts nicely for the time evolution 
of this quantity for all system sizes for t/L < 10 4 . This 
is a second puzzle we will try to address further in the 
following sections as well as in [ffsl . 



3. Mean Hamming distance 

We also calculate the Hamming distance of two succes- 
sive shear bands which is defined as the number of sites 
(or bonds) by which two consecutive shear bands differ. 
We denote this distance by d. The natural normalization 
is to divide this distance by the total length of the path, 
L x . As we will see this quantity is very useful in charac- 
terizing the time evolution of the localization process. 




JO" L-J .1 . . .1 . . .1 . .1 . .1 .1 .1 I L 

10 " 3 10 1 10 1 10 3 10 5 

t/L 

FIG. 5: The difference of the mean density of the shear band 
from its asymptotic value 0.5 vs. time. Notation and system 
sizes are the same as in Fig. 0. 
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t/L 

FIG. 6: The average Hamming distance versus time. The 
same system sizes were scaled together as on Fig. ^. The an- 
alytical prediction l/(t+l) is plotted over the data. Note 
that scaling with system size displays systematic corrections 
for larger systems. 



Fig. |6| shows that the mean Hamming distance is close 
to unity (i.e. two consecutive paths do not overlap at all) 
at early times, and decrease towards for t/L ^> 1. We 
recall that when the distance is equal to 0, then the two 
consecutive conformations of the shear band are iden- 
tical, in spite of the total renewal of random densities 
along them. This indicates that the shear bands have 
a tendency to remain more and more persistent as the 
system "ages" . We will analyze further the complete sta- 
tistical distribution of the Hamming distance later in this 
paper. 
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7: Numerical Estimation of the cumulative shearing 
a cum- Figure a) represents the number of times (n a ) a site 
y was active up to time t = 1000 in a cross-section of a 128 
by 128 sample. Figure b) is the cumulative representation of 
a): ^,j =0 n a (j)- The dashed line indicates the homogeneous 
case. 
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FIG. 8: The cumulative shear corrected by the average dis- 
placement for a system of size 128 by 128. As one can observe, 
this quantity encodes the history of the process; a dip in the 
profile indicates the presence of the shear band and the depth 
of the dip is indicative of the amount of time it has spent in 
that location. For example, in the beginning when successive 
shear bands are distinct, every site is visited approximately 
equally and the profile has no deep peaks or valleys. After 
this, at t = 5000, the first shear band gets localized at around 
x = 100. This lasts until about t = 20, 000, and then it jumps 
to x = 40. After spending some time here it jumps back, close 
to its previous position (2 10 4 < t < 5 10 4 ). 



An experimentally relevant quantity is the cumulative 
shearing denoted by a cum . The numerical procedure we 
follow to obtain this quantity in our model is the following 
(see Fig. 0): We mark a line in the y-direction (see Fig. 
0). We measure the total activity n a {y) along the line, 
i.e. the number of instances when the shear band went 
through a point y on the line. From this, we define 

y 

Vcum(y) =y ^n a [j). (2) 

3=0 

By definition a cum {L y ) = t, since at every instant, the 
shear band has necessarily to pass through one or the 
other site on a cut along the y axis. The fluctuations 
of cr cum (y) about its mean value (t/L y )y then reflect the 
inhomogeneity of the shear process within the sample. 
In Fig. we track the time evolution of cr cum (y), after 



subtracting out the mean value. As indicated, a snap- 
shot of this quantity encodes the history of the process 
of shearing in this system. 



B. Early time regime 

It is apparent from Figures ||, || and ^ that the ini- 
tial behavior of the model is very different from the late 
stages. In the former regime, the average distance d (Fig. 
^) is equal to the system size indicating that successive 
shear bands do not overlap at all. In other words there 
is an effective strong repulsive interaction between them. 

The density on the shear band (Fig. ||) provides an 
explanation for this behavior. From the directed poly- 
mer [fill picture it is known that the first shear band 
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has a mean density of gsB(t — 0) = e* w 0.22 on the 
tilted square lattice. The state of the shear band after 
the densities of the sites along it have been changed, is 
independent of its previous state and the mean density 
is now 0.5 (since it is just an average of L x indepen- 
dent random numbers taken from the uniform distribu- 
tion [0,1]). This implies that at the next instant, the 
chosen shear band will have a very low probability of 
sharing bonds/sites with the previous one, since there 
will still be many spanning paths with density smaller 
than 0.5. Thus, until all the sites are visited at least once, 
successive shear bands have few sites/bonds in common 
(d = L x in Figs. ||). Also in this regime, since the path 
sweeps all sites through the cumulative shearing appears 
to be homogeneous (Fig. || first plot). As mentioned in 
the introduction, this uniform shear strain is consistent 
with the experimental observations that no well-defined 
shear-bands persist over observable time scales in loose 
samples. 

The characteristic time to build up correlations is set 
by the sweeping through the sample, i.e., it is given by 
L y . This is the reason why in the early time regime the 
plots in Figs. 4-6 can be scaled together with L y . 



C. Localization 

The above described behavior is drastically changed as 
time goes on and the shear band gets localized for very 
long times at the same location. 

As the average density increases, the probability of 
choosing a minimal path with density less than 0.5 starts 
decreasing, and thus for the new path, it becomes more 
favorable to overlap to a greater and greater extent with 
the previous one. As a result, the density of the region 
where the path is located, is further increased, hence 
trapping the minimal path in canyon-like structures sur- 
rounded by extremely high density regions (see Fig. ^ 
b-d). 

The early repulsive interaction is thus now inverted 
to an attractive one as is shown by the rapid 1/t like 
decrease of the Hamming distance (Fig. ||), as well as 
the decrease of 0.5 — qsb and of (g) from its initial value 
(Figs. H and [I| respectively) in this regime. 

There are a number of consequences of the localiza- 
tion. First, as time goes on, the shear band (which was 
not visible at all in the density map of the system (see 
Fig. |3| a)) becomes more and more apparent until fi- 
nally it becomes localized at a given position for macro- 
scopic times. This is in accordance with experiments and 
with the critical state concept ||, E(J. Secondly, since 
the same path for the shear band is chosen most of the 
time, its density saturates to its asymptotic value 0.5 
(Fig. ||). Again this is consistent with the experimen- 
tal observation that the density within the shear band 
tends to achieve a well defined value somewhat smaller 
that the rest of the medium for dense granular media. 
Simultaneously, the shear stress saturates to a constant 



value since this is imposed by the shear band itself [Ej. 
However, we note that in our model the global density 
of the system continues to increase in time, albeit very 
slowly. This extremely slow trend may well be out of 
reach experimentally. However, it is known jj] that the 
shear stress saturates much faster than the volumetric 
strain (the average density in our case) which is clearly 
justified by the numerical results. 

On the cumulative shear which is a straight line with 
small statistical fluctuations in the early time regime, 
there appear step-like structures, indicators of progres- 
sively more persistent localization (Fig. ||). However, 
this localization is not everlasting since the shear band 
may perform big jumps to other local minima. This can 
be seen on the series of cumulative shear curves in the 
form of certain steps disappearing and others becoming 
more prominent. This prediction could easily be tested 
experimentally. 



D. Systems with different aspect ratios 

All changes take place along the shear band, which is 
aligned along the x direction, and thus we may antici- 
pate that the x and y direction will play different roles. 
Therefore in this section, we study the influence of the 
width and length of the system. In what follows, we use 
the terms long for samples with L y < L x and wide in the 
opposite case (L y > L x ). 

With very long and wide systems we are able to sepa - 



rate the two kinds of dynamics described in Section IV C . 
If one considers a long sample [L v /L x is small), we expect 
the time evolution to be independent of L x for all quanti- 
ties of interest since the lattice can be split into subparts 
placed in series. Thus one may expect the large jumps 
to disappear and the average density to scale solely with 
L y (Fig. ||). A wide system, on the other hand, might be 
expected a priori to behave like a number of competing 
subsystems. 



1. Long systems 

On Fig. ^| we have plotted the time dependence of the 
density in long samples with L y = 4 (lower curves) and 
L y = 8 (upper curves). The t/L y scaling is excellent in 
both cases. However, the densification obeys a different 
time evolution for different L y . The rate at which the 
density increases is slower as the width increases. 



2. Wide systems: Breakdown of ergodicity 

Wide samples can be considered as subsystems placed 
next to each other and coupled in parallel. In contrast 
with the previous case, we will see that the evolution of 
the different subsystems cannot be accounted for by a 
simple average. 
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FIG. 9: Long samples with width L y — 4 (lower curves) and 
L y = 8 (upper curves). Three different lengths were used in 
both cases L x /L y = 2.5, 5 and 10. 
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FIG. 10: Time dependence of the density difference from its 
asymptotic value for wide samples. The length of the system 
along the shear band is L x = 5 the widths are L y — 5, 10, 20, 
40, 80 and 160 from bottom to top, respectively. 



sider an extreme version of such a breakdown of ergod- 
icity. During a first stage, up to t/L y of order 1, the 
activity is evenly spread over the system. Then we as- 
sume that after such a time, the activity remains confined 
in a subsystem of size L x x I. All other subsystems are 
assumed not to be visited by the shear band, and thus 
their density is quenched at their value reached at the 



onset of localization, go, at time to 
density will thus obey 



6L y . The global 



{g){t) _ (1 - Qo)(L y -i) + fi(t - 6(L y - £))£ 



where ft(t) = 1 — (g)(t) describes the densification of 
the representative cell of size L x x £ where the activity is 
confined, and thus 1 — go = ft{6l). 

In this crude scenario, we note that the global den- 
sity does not converge to 1 as time goes to infinity, but 
rather remains stuck at a value such that (1 — (g)(t)) — » 
fi(6£)(l —t/Ly). In more quantitative terms, we tried to 
carry out such a procedure, and indeed for a fixed L Xl it 
is possible to account for the time evolution of systems of 
different width using go as free parameter and calculating 
I from the L y dependence. It turns out that I changes for 
small values of L x but becomes constant {£ ~ 30) above 
the system size of L x ~ 30. The test of this analysis can 
be seen on Fig. |ll|. However, the asymptotic density 
turns out to depend on L x . 

The conclusion is that although, such an extreme mod- 
eling of the localization is able to capture part of the 
strong size effects observed numerically, it is too crude to 
provide a quantitative account of the densification. The 
hierarchical lattice provides us with a convenient case 
where an analytical investigation of this breakdown of 
ergodicity can be made. It is shown in Ref. p8| , that 
the local "age" distribution assumes a multifractal distr- 
bution whose spectrum can be computed exactly. This 
property can then be used to provide an expression of 
the density evolution in time. 



In the case of the wide systems the same plot as before 
(Fig. [h] a) shows no t/L y scaling for small system sizes. 

We could imagine the following construction: suppose 
we split a given wide system into two subsystems of size 
L x x (Ly/2). Provided L y is large enough, we can ignore 
the interaction between the two subsystems, and thus 
we could study independently the time evolution of both 
subsystems. Now, if we merge them again, we realize 
that the only reason why the resulting densification could 
differ from the measurement on the separate subsystems 
is that the time t\ the shear band has stayed in subpart 
1 is far from being equal to t/2. In other words, the 
breakdown of the data collapse of 1 — (g) vs. t/L y is a 
breakdown of ergodicity. This is naturally associated to 
what we termed "localization" earlier. 

In order to make this concept more explicit, let us con- 



E. Time evolution of the Hamming distance 
distribution 

We study here the distribution of Hamming distances 
as a function of time, P(d, t) for both lattices. This quan- 
tity is the analogue of an "avalanche distribution" , such 
as is usually studied in self-organized critical(SOC) sys- 
tems. As we will see below, this quantity does indeed de- 
cay as a power-law like in SOC systems. However, since 
a steady state is never reached, the power-law decay is 
multiplied by a time-dependent prefactor. 

The distribution P(d, t) is shown in Fig. [l2| At early 
times, this quantity is peaked around the maximum value 
(d = L x ) while in the localized regime, it becomes peaked 
at the minimum value [d = 0). This corresponds to the 
transition from repulsive to attractive effective interac- 
tion between consecutive path conformations. 
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FIG. 11: The same data as in Fig. [jlj with a different, effective 
density g* = Qo(L y — I). Note that the plot displays nice 
scaling as well as a clean power law decay over at least three 
decades. However if the simulation is continued further the 
average density increases over g* . 



At fixed (large) times, the distribution P(d,t) decays 
as a power-law of d, as can be seen in Figure ftq. The 
measured exponent is 



P(d,t) oc dT 



(4) 



in addition to which there exists a peak at d = the am- 
plitude of which varies significantly with time. The decay 
of average dasl/t (see Fig. ||) found earlier, implies that 
the time dependence of the d ^ part is p(d, t) oc 1 / (£<i 3 ) . 
Thus, including the different scalings with L x and L y , we 
obtain finally the asymptotic form 



(5) 



Here again, the hierarchical lattice allows us to com- 
pute this distribution analytically (for L x — L„), and we 
find that a similar behavior is obtained (Lai. 



F. Roughness exponent of the shear band 

From the directed polymer analogy we know that the 
shape of the shear band is self-affine with an exponent of 
C = 2/3 for infinitely large systems. It is an interesting 
question whether this property of self-affinity is conserved 
in the time evolution of our system. We have investigated 
this question and have found self-affine scaling for all 
times, albeit with a time dependent Hurst exponent (Fig. 

We have estimated the value of C by measuring the 
width of the shear band, wi,(t), for different system sizes, 
L x L. The width is defined as the standard deviation of 
the y coordinate of the active path. The latter is expected 
to scale as wz,(t) oc L> for a self-affine object. To estimate 




FIG. 12: The distribution of the Hamming distance on a 
64x64 square system for different times. The dashed lines 
indicate the mean (d). 
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FIG. 13: Scaling plot for the distribution of the Hamming 
distance P(d, t) vs. d. The three curves are for three square 
system of size: L = 64 at t = 10000, L = 128 at t = 1000 
and L = 256 at t — 500. The straight line has slope —3 
indicating that the decay of P(d, t) with d is a power-law. 
Jumps of order of half the system-size however, seem to have 
an enhanced probability. 
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FIG. 14: Estimated self-affine exponent as a function of time 
obtained through Eq. ^| for system sizes varying from L — 64 
to L = 512 . 

C, the roughness exponent, we have computed the ratio 
of two such widths for lattice sizes differing by a factor 
of 2, and used the following estimate 

= log {w 2L (t)/w L (t)) 
<■ log(2) 1 > 

where (t) is the width of a shear band inaLxL lattice 
at time t. The results obtained for the tilted square lat- 
tice can be seen on Fig. [l4[ It starts from £(i = 0) = 2/3 
as expected from the directed polymer result and has an 
asymptotic value of £ ~ 0.8. 
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FIG. 15: Time dependence of the difference of the average 
density from its asymptotic value for starting densities with 
initial density ranges [Qi n u : 1] from top to bottom respec- 
tively: Q ini t = 0.2, 0.3, 0.4, 0.45, 0.5, 0.55, 0.6, 0.7, 0.8. The 
system size is 64. 



expectation value of the refreshing density distribution 
= 0.5. After the first regime, as all small values are elim- 
inated, the system effaces the initial condition almost 
entirely. 

On the other hand, starting from initial conditions 
with Qinu > 0.5 we largely eliminate the possibility of 
big jumps. The shear band fluctuations are now very 
small, involving changes in a very few sites, and hence, 
the density change of the sample is extremely slow. 



G. Systems with high initial densities; relevance of 
initial conditions 

We have so far studied the situation when a very loose 
granular medium compactifies under shear, while simul- 
taneously a shear band gets quasi-localized in the system. 
It is also of interest to study samples with higher initial 
density where it is known experimentally that a shear 
band is localized from the very beginning. In this section 
we study the interesting crossover to that state from the 
previously described dynamics. 

We choose initial densities from the interval [Qinit ■ 1] 
with a uniform distribution and varying . In Fig. [l5| 
the time evolution of the difference of the average density 
from its asymptotic value is plotted using different initial 
conditions. Since all previous arguments hold we assume 
that in finite systems the asymptotic value of the average 
density is 1 irrespective of the initial density. 

The striking result of these simulations is that if Qinit < 
0.5, all curves coincide in the decreasing regime. How- 
ever, if Qinit > 0.5 a different time evolution is observed 
for large times. Thus we can assume that in the early 
time regime, when the shear band is swapping uncorre- 
latedly, it visits all sites that have a value less than the 



H. Summary of the numerical results 

We have seen that the system densifies with time so as 
to approach a unit density, i.e. the maximum available 
density from the distribution used to refresh the sites. 
The kinetics of the densification is slow (slower than any 
power-law). Moreover, after a first transient where the 
reduced time t/L accounts for the L dependence, the 
compaction process depends on the system size in a non- 
trivial way. The width of the system is the parameter 
which really controls this anomalous behavior, signaling 
that the competition between parallel paths may some- 
how play a key role in this breakdown of system size 
rescaling. This competition is a subtle one however, ly- 
ing somewhere in between complete localization of the 
path (which, as we saw in section IV B, is too crude to 
mirror the actual scenario) and complete derealization 
(which, as mentioned earlier, accounts only for the early 
time behaviour). 

The density maps display an interesting organization 
of "canyon-like" paths with density much lower than their 
immediate surroundings (where the density approaches 1 
quite uniformly). Moreover, large regions are left quies- 
cent, being systematically avoided by the minimal paths. 
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This contrast of high and low activity within the same 
system is at the heart of the breakdown of ergodicity ob- 
served after an initial transient. As remarked however, 
the distribution of Hamming distances between consecu- 
tive paths displays a somewhat simpler behavior where 
the role of the width and the length of the system can be 
simply accounted for. 

V. CONCLUSION AND DISCUSSION 

Though very simply defined, our model seems to cap- 
ture some essential features of granular shear and pro- 
vides at the same time several predictions. The model 
demonstrates the self-organized mechanism of the local- 
ization of the shear band in loose granular materials. As 
the sample ages, very high fluctuations in density ap- 
pear where we can observe some kind of screening ef- 
fect: more resistant regions of higher density, protect the 
looser ones. This model also gives an insight into a dy- 
namics that exhibits very non-trivial system size effects. 

Our stochastic model makes an attempt to describe 
the large strain behavior of sheared loose granular mat- 
ter on a mesoscopic level. The rules of the model do not 
include any dependence on the total amount of shear im- 
posed on the medium, nevertheless, a constant friction 
angle and slow densification is observed — a property 
referred to as "ageing" — which reproduces the exper- 
imental results qualitatively jy]. By construction, the 
strain takes place through local shear bands which ini- 



tially travel throughout the medium homogeneously (and 
hence produce a uniform shear), but which progressively 
become more permanent giving rise to more steady shear 
bands, a feature also observed experimentally p3| . Our 
model reproduces further features seen in experiments 
and numerical simulations, including the high frequency 
fluctuations of the local shear |(| . 

In addition, we predict a complex self-organization of 
these shear bands, displayed in the inhomogeneities in the 
local density. This feature can be studied experimentally, 
in particular through the use of X-ray tomography, to 
access the local density of a sheared medium. The use 
of tracer particles could also be helpful in identifying the 
inhomogeneous ageing and localization of the shear bands 
as well as their sudden changes. 

Most of the results presented here for the Euclidean 
lattice are closely mirrored by the results on the hierar- 
chical lattice studied in |18|j . The recursive topology of 
this lattice allows a quantitative analytical understand- 
ing of many of the quantities studied numerically in this 
paper. This includes elucidating the mechanism for the 
breakdown of ergodicity, the slow density evolution, as 
well as the behaviour of the of the Hamming distance at 
late times. 
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